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Abstract 

We consider the problem of recovering a set of correlated signals {e.g., images from different viewpoints) 
from a few linear measurements per signal. We assume that each sensor in a network acquires a compressed 
signal in the form of linear measurements and sends it to a joint decoder for reconstruction. We propose 
a novel joint reconstruction algorithm that exploits correlation among underlying signals. Our correlation 
model considers geometrical transformations between the supports of the different signals. The proposed 
joint decoder estimates the correlation and reconstructs the signals using a simple thresholding algorithm. We 
give both theoretical and experimental evidence to show that our method largely outperforms independent 
decoding in terms of support recovery and reconstruction quality. 



1 Introduction 

The growing number of distributed systems in recent years has led to an important body of work on the efficient 
representation of signals captured by multiple sensors. Recently, ideas based on Compressed Sensing (CS) [Tl[2] 
have been applied to distributed reconstruction problems f3| in order to recover signals from a few measurements 
per sensor. When signals are correlated, a joint decoder that properly exploits the inter-sensor dependencies 
is expected to outperform independent decoding in terms of reconstruction quality. Very often, the correlation 
model restricts the unknown signals to share a common support. Using this correlation model, the authors in 
and [4] propose decoding algorithms and show analytically that joint reconstruction outperforms independent 
reconstructions. In many applications, this correlation model is however too restrictive. For example, in the case 
of a network of neighbouring cameras capturing one scene or seismic signals captured via different sismometers, 
the supports of the signals are quite different even if components are linked by simple transformations. 

In this paper, we adopt a more general correlation model and build a joint decoder that recovers the 
unknown signals from a few measurements per sensor. We assume that the unknown signals are sparse in a 
redundant dictionary V, and not necessarily in an orthonormal basis [5l[6]. We denote the components of the 
dictionary as atoms. We assume that the support of each view j is related to the support of a reference view 
by a transformation T* . The transformation T* could be for example a translation function. Using the given 
correlation model, we build a joint decoder based on the thresholding algorithm ,5j and prove theoretically that 
it outperforms an independent decoding method in terms of recovery rate. Moreover, we show experimentally 
that the proposed algorihm leads to better reconstruction quality. 



2 Problem formulation 

We consider a sensor network of J nodes. Each sensor j acquires M linear measurements of the unknown signal 
yj G M.^ {M < N) and sends it to a central decoder . The role of the decoder is to estimate the unknown signals 
y = {yj}j^i. By denoting S — {sj}j^^ the set of compressed signals acquired by the sensors and A — {Aj}^^^ 
the sensing matrices, we have: 

(1) 

Mxl MxN Nxl 

In the rest of this paper, we use independent sensing matrices with Gaussian i.i.d entries. Specifically, \/M{Aj)m,n 
follows a standard Gaussian distribution, for any m,n,j. 

We assume that the unknown signals yj G are sparse in some dictionary V that consists of K atoms 
and denote by $ = [ipi \ . . . I^pk] its matrix representation. Formally, we have yj = ^Cj, where Cj is a vector of 
length K with at most S non zero components and S < N. By denoting the support of yj with A* (i.e., the 
set of S atoms corresponding to the non zero entries of Cj), yj can be written as follows: 

^ = (2) 

Nx-l NxS Sx.1 

where is the restriction of to A* and Xj corresponds to the non zero entries of Cj. 



We adopt the following correlation model: for any j e {!,..., J}, the set of atoms in A* can be obtained 
from A^ by applying a transformation T* : V ^ V. This can be written as: T*{/S.i) — A* (we consider that 
is the identity). In our problem, the vector of transformations T* — {T*}'j^^ is unknown. However, we assume 
that we are given a finite set T of candidate transformations vectors and that the correct vector T* belongs to 

r. 

Considering the above correlation model, we address the following problem: Given the compressed signals 
iS, the sensing matrices the sparsity 5, the dictionary $, and the set of candidate transformations vectors 
T, estimate the unknown signals y [i.e., supports {A*}^^]^ and coefficients {xj^j^^) using a small number of 
measurements per sensor M . 



3 Joint thresholding algorithm 

We propose a solution to the problem formulated in the previous section. Our proposed decoder extends 
the simple thresholding algorithm [5 to multiple signals. This choice is motivated by the low complexity of 
thresholding algorithm with respect to other decoding methods [7]. Our joint decoder represents an efficient 
alternative when the signals are simple {i.e., they have very sparse representations in the dictionary) and the 
number of sensors is fairly large, so that other decoding methods become computationally intractable. 

The Joint Thresholding (JT) decoder exploits the information diversity brought by the different signals 
to reduce the number of measurements per sensor required for accurate signals reconstruction. It groups 
the measurements obtained from each individual signal and precisely estimates the unknowns (A*,T*) (or 
equivalently all the supports {A*}^^;^). 

JT obtains an estimate (Ai,r) of (A^,T*) by maximizing the following objective function, which is called 
the score function: 

J 

vI/,(Ai,r) = ^ ^ s,-A,^, (3) 
j=i veT,(Ai) 

where Ai and T denote respectively the support of the reference signal and the vector of transformations 
variables and the operator • denotes the canonical inner product. The use of ^'^ as the objective function is 
justified by: 

(A;,f) =argmax4',(Ai,r) 

(Ai,T) 

« argmax*2^(Ai,r) 

(Ai,T) 

= (AI,T*), 

where 4',(Ai, T) = ^.ti E^eT,(A, )Vj-V = (E denotes the expected value of a random variable). Indeed, 
if both assumptions given in Eq.Q and Eq.([5| holcQ ^'y(Ai,T) is maximal for Ai = AJ and T = T* . Besides, 



for large values of M, 5's(Ai, T) concentrates around its average value 5'j,(Ai, T) (Lemma 4.1 ). The description 

of JT algorithm is given in Algorithm [T] 

In words, the JT algorithm calculates for each transformation vector T e T the vector dx, whose entries are 
J 

given by drli] = *i ' ■^j'^jiVi)^ for ^ < i < K . Then, the largest S elements in are summed and assigned 

to ^'s(Ai,r). Estimated quantities {Ai,T} are updated if 5's(Ai,T) achieves a higher score. Knowing the set 
of supports, we deduce coefficients Xj by computing the least squares solution to equation Aj<^^Xj = Sj. 



4 Theoretical analysis 

Our theoretical analysis focuses on the performance of JT in finding the correct supports. Hence, we will not 
address the quality of the estimated coefficients {xj}^^^. In particular, we focus on the analysis of the recovery 
rate R defined as the total number of correctly recovered atoms (in all signals combined) divided by the total 
number of atoms (i.e., SJ). 
We assume the following: 



• For all j e { 1 , . . . , J} , A* can be recovered entirely by applying the thresholding algorithm on Uj . Formally, 
there exists 77 > verifying: 



inf 



■ f 



> sup 

VSA* 



■ f 



V, 



(4) 



^The assumptions are discussed in the next section. 
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Algorithm 1 Joint Thresliolding (JT) Algoritlim 



Input: compressed signals {sj}, sensing matrices {Aj}, sparsity S, dictionary candidate vectors of trans- 
formations T 

Output: estimated signals {yj}, support Ai and vector of transformations T 

1. Initialization: (A^,?, $) ^ {0,0, -oo). 

2. For every T e T 

2.1 Build the vector dx of length K in the following way: 

.7 

dT = ^(A,r,.($))%-, 

whereT,i^)^[Tj{cp,)...T,{cpK)]. 

2.2 Keep the largest S entries in dx and set the other entries to zero. The positions of the non zero entries 
in dx give the indices of the estimated support Ai of the first signal. 

2.3 Calculate the score ^'s(Ai,r) by summing the S non zero entries of dy. 

2.4 If ^'^(Ai,r) exceeds $: update (A!,f,$) ^ (Ai, T, ^',.(Ai, T)) 

3. Build the coefficients vector for all j £ {1, . . . , J}: 

where (•)+ denotes the pseudo-inverse operator. Note that Aj is obtained using the correlation model: 
a; =T;(A;), for j > 2. 

4. Obtain signals estimates: 



where A* is the complement of A* in V. As this condition is practically hard to verify, a sufficient 
condition involving the coherence of the dictionary is given in -5, Eq.(3.2)]. 



• All the atoms in the supports have positive inner products with the corresponding signal: 



(5) 



The assumption in Eq.Q is reasonable since we cannot hope to recover the supports using JT unless the 
thresholding algorithm correctly recovers the supports when applied on the full signals yj . Assumption in Eq. ([s]) 
is a technical one and is used in the proof of our main theorem. Intuitively, it guarantees that ^'y(Ai,r) = 
E^'s(Ai,r) is maximal when (Ai,T) = {Al,T*). This assumption can be achieved by adding the inverse of 
the atom in the dictionary ((/? — ^ —(p) when the inner product is negative. 

The main ingredient we will use in our analysis is the concentration of j ^j^j '^j'^^j ai'ound its average 



value jJ2j=i' 
lemma: 



for any set of vectors {uj}j^-^ and {vj}j^^ of length A^. This is shown in the following 



Lemma 4.1. Let {uj}i<j<j and {vj}i<j<j with Uj,Vj £ E , such that \\uj\\2 < and \\vj\\2 < -Bu for all 
j G {1, . . . , J}. Assume that {>lj}i<j<j are independent random matrices of dimension M x N , with iid entries 
following Af{Q, jj). Then, for all t > 0, 



^AjUj-AjVj 



> T < 2 exp 



C^BIBI + C^tB^BJ ■ 



(6) 



with Ci 



and C2 = 2%/2e 



The proof of this lemma can be found in Appendix A. 

Our main theoretical result is given in the following theorem: 

Theorem 4.1 (Recovery rate of JT). Let R be the recovery rate of JT defined by: 



E"^i|A*nA, 

SJ 
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Then, for any < a < 1; 



-CMJrfa^^] , (7) 



where rUy = minj \\y.j\\2, My = maxj ||j/-,||2, C = y^i; + 4e\/2j . 
The proof of this theorem can be found in Appendix B. 

For simphcity, we consider the common case where aU the signals have the same energy [my = My). Theorem 
4.1 shows that for sufficiently high values of J, the recovery rate is mainly governed by MJ, ry and |T|. The 
dependence on MJ (i.e., total number of measurements) follows our intuition as JT combines the measurements 
of the different sensors to perform the joint decoding. Increasing the total number of measurements leads to 
a better recovery rate. The quantity rj hides the dependence of R on the signal characteristics and model. 
For clarification, the following inequality provides a lower bound on 77, in terms of sparsity, coherence of the 
dictionary and ratio between the lowest to largest coefficients: 

(telu -Mi(^-1)-Mi(^))' 



if' > min ■ 



3 Sil + fiiiS-1)) 

where ni defines the cumulative coherence (Babel function) as defined in |8] and |xniin,j| is the absolute value 
of the smallest coefficient in vector Xj. Note that if $ is an ONB, /ii = 0. The proof of this inequality is very 
similar to the proof of Corollary 3.3 in [S]. 

Another key parameter is the number of candidate vectors of transformations |T| which grows with J. In 
the following corollary, we provide a lower bound on the number of measurements needed per sensor to reach 
asymptotically a perfect recovery rate in the following two cases: (1) T grows slowly with J, (2) T grows 
exponentially with J. 

Corrolary 4.1 (Asymptotic behaviour of R). Let < a < 1. 

1. If\T\ is a subexponential function of J , then, as long as M > 1, V{R > 1 — a) converges to 1 as J ^ +00. 

2. If there exists /3 > such that \T\ ~ e^'^ , then, as long as M > q^i^i ^[R > 1 — a) converges to 1 
as J ^ +00. 



Proof. The proof is a direct application of Theorem |4.1| □ 
The growth of |T| is related to the degree of uncertainty on the correct transformation vector T* . For 



example, if T* is known in advance, then |T| = 1 and Corrolary 4.1 guarantees an arbitrary high recovery rate 
with only one measurement per sensor when J +00. This result remains valid as long as |T| <C e^"^ for all 
/3 > 0. However, if T* is completely unknown and transforms between pairs of signals are independent, |T| grows 
exponentially with J and we will need more measurements per sensor in order to recover the correct support 
estimates (consider the example where (a) the number of candidate transformations between each sensor j > 2 
and the reference signal is equal to / ; (b) Tj is independent of Tj_i, then: |T| — l''~^). 

Unlike independent thresholding which has a constant recovery rate in function of J, previous results show 
that the recovery rate of JT increases by augmenting the number of sensors J. Thus, in large networks, JT 
requires less measurements per sensor than independent thresholding for a fixed target recovery rate provided 
that |T| has a controlled growth. 



5 Experimental results 
5.1 Greedy JT 

The JT algorithm, as described in Section [3j performs the search over all candidate transforms in T. This 
can be very costly in terms of the computational efficiency, especially for a large number of correlated signals. 
Thus, instead of performing a full search, we greedily look for the relevant transformations. The Greedy Joint 
Thresholding (GJT) algorithm is described in Algorithm [2] 

Even though this algorithm has a lower complexity than JT, the price to pay is a less robust transformation 
estimation process: in the early stages of the algorithm (y <C J), the selection of the transform is based on a 
small number of signals V . If in addition the value of M is small, this may lead to uncorrect estimation of the 
transformations and thus wrong support estimates. 

In Figjl] we plot the percentage of uncorrect estimated transforms with JT and Greedy JT in function of 
M, for a randomly generated image. For M > 80, the performance loss with Greedy JT is relatively small with 
respect to the gain in complexity. Thus the penalty of using the greedy algorithm is small in practice. In the 
following, we examine the performance of Greedy JT on synthetic images and seismic signals. 
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Algorithm 2 Greedy Joint Thresholding (GJT) Algorithm 



Input: compressed signals {sj}, sensing matrices {Aj}, sparsity S, dictionary candidate vectors of trans- 
formations T 

Output: estimated signals {yj}, support Ai and vector of transformations T 

1. Initialization: T ^ 1 (identity). 

2. For every F e {2, . . . , J} 

2.1 Set the values (Ai,7V,$) ^ (0,0, -oo). 

2.2 Let TV denote the possible transformations between signal 1 and signal V. 

2.3 For each element Ty G TV do 

2.3.1 Let T 4- [f ,Ty] (i.e., Ty appended to f) 

2.3.2 Compute: 

V 

2.3.3 Keep the largest S entries (set the other entries to zero). The positions of the non zero entries in dx 
give the indices of the estimated support Ai. 

2.3.4 Calculate the score ^'s(Ai,r) by summing the S non zero entries of dx 

2.3.5 If ^'^(Ai,r) exceeds $: update (A!,IV,$) ^ (Ai, Ty, *,(Ai, T)) 

2.4 Update the estimate of the vector of transformations: T [T. Ty]- 

3. Perform steps 3. and 4. in Algorithm [l] 




5.2 Synthetic images 

We construct a parametric dictionary where a generating function undergoes rotation, scaling and translation 
operations to generate the different atoms in the dictionary We use the Gaussian g{x,y) ~ e^^ as the 
generating function. The atoms in the dictionary are characterized by the rotation angle 9, scales Sx and Sy 
and translations and ty. If {X,Y) denotes the transformed coordinate system: 

^ _ ix-tx)cos9- {y-ty)sm9 

Sx 

_ {y -ty)cos9 + {x~tx)sin9 

Sy 

the atom gp with parameters p — {9 , Sx, Sy,tx,ty) is given by: 

gp{x,y) = pg(A,r), 

where p is the normalization constant. 

The dictionary is generated for images of size iV = 32 x 32 = 1024, with the following parameters: 9 G 
[0 : ^ : 7r],Sa: = {2,4},Sj, = {1/2, 1}. Every atom is shifted in pixels of odd coordinates, so the full dictionary 
contains 6144 atoms. 

The support of the reference image and coefficients are chosen in order to verify the conditions in Eq.Q and 
Eq.(l5]). The remaining images have been obtained by applying global translations on the atoms of the reference 
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image, under the constraint that all atoms in images belong to V. We assume that the transformations are 
independent from one another and that there are 9 candidate transformations for any image. Thus, |T| — O'^^^. 
Figj2]illustrates the recovery rate and MSB of Greedy JT and independent thresholding for a randomly generated 
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Figure 2. Recovery rate and Mean Squared Error of Greedy JT (GJT) and Independent thresh- 
olding (IT) in function of J. Simulation setup: 10 independent trials, S = 5, M = 150, 
A'^ = 1024, Gaussian sensing matrices. 

image. Recovery rate is defined in Section |4j For a given J, the calculated MSE represents the averaged MSE 
calculated over signals {!,..., J}. We see that Greedy JT outperforms independent thresholding in terms of 
recovery rate and image quality, especially for high values of J ( J > 20). Thus, although |T| grows rapidly with 
J, our joint decoding approach is significantly better in practice in terms of support recovery. 



5.3 ID seismic signals 

Seismic signals captured at neighbouring locations typically follow the correlation model proposed in this paper. 
Figjs] (a) , (b) represent two seismic signals that are obviously correlated as the second signal is approximately 
a shifted version toward the front of the first signal. We use the following sparsifying dictionary, which consists 
of Gaussians modulated with sinusoids: 



9(t,s,ui){x) = pexp 



where p is the normalizing constant. The translations t are chosen uniformly from 1 to iV with step size 10 such 
that the coherence of the dictionary is not too high. Scales s take values in {4, 8, 16} and u! varies from 2 to 10 
with step 2. For each set of parameters (t, s,aj), g(t,s.uj) ^-^id —g(t.s,uj) ^-re included in the dictionary. Figjs] (c) 
and (d) illustrate the estimations of signal number 2 obtained with only 15% of the measurements respectively 
using independent thresholding and JT algorithm. Note that as J = 2 in this example, Greedy JT and JT are 
equivalent. Visual inspection and calculated MSEs confirm the superiority of joint decoding using JT algorithm 
over independent thresholding in terms of reconstruction quality. This experiment shows that joint decoding 
using JT provides significantly better quality signals even when the number of correlated signals is low (J — 2). 



6 Conclusion 

In this paper, we have proposed an efficient approach for the joint recovery of correlated signals that have been 
compressed independently. Our solution is novel with respect to the state of the art work due to the particular 
geometrical correlation model based on the transformations of the sparse signal components. Mathemati- 
cal analysis and experimental results demonstrate the superiority of our recovery algorithm over independent 
thresholding. JT is namely applicable for decoding simple multiview images, seismic signals or any other set 
of correlated signals satisfying the geometric correlation model. A promising future direction is to use JT for 
correlation estimation along with a more sophisticated recovery algorithm for the reconstruction. 
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Figure 3. Seismic signals (a) j/i and (b) j/2 captured at two neighbouring locations. Estimation 
of j/2 using (c) independent thresholding and (d) JT. Simulation setup: J = 2, N = 1000, M = 
150, S = 50, |T| = 3, Gaussian sensing matrices. This experiment was conducted 200 times and 
we obtained MSEu = 0.0031 and MSEjj = 0.0025. 
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Appendix A: Proof of Lemma |4.1 



This proof is inspired from the proof of [SJ Lerama 3.1]. 

Let {gj)m,n be a random variable foUowing the standard gaussian distribution such that {Aj)„i^n — ^7^^ {gj)r 
We have, for any set of vectors Uj,Vj in M^: 

^ J ^ J M N N 



J 



j=l m=l k=l 1=1 



N N 



Let Y,,™ = X) Y.{9j)m,k{9])m,i{uj)k{vj)i- As (gjOm.fc is independent from igj)m.h we have EY^_,„ = uj-Vj. Let 

fc=l /=! 

Zj,m = ~ E^j.m = J2 {9j)ni,k{9j)m.i{uj)k{vj)i + J2iii9j)ni,kf " l){uj)k{vj)k- By definition, for any j and 

fe^^/ A; 

m, .^j.m is a Gaussian chaos of order 2. 

Observe that the probability we wish to bound can be expressed in terms of Zj^m- 



jVj-Uj - v. J 



J M 



M 

j — 1 7n—l 



> T 



J M 




E E ^j>" 




j — l m— 1 





As Zj „i is a Gaussian chaos of order 2, Bernstein's inequality on the sum of zero mean independent random 
variables with a certain moment growth is applicable (notice that the independence assumption is satisfied in 
our case as entries are iid and {Aj}'j^^ are independent). For more details about the theorem refer to Theorem 
A.l in [5]. We get: 



.7 M 

E E ^j^™ 
j=i "1=1 



>tMJ\ < 2 exp 



1 JMt2 

2 W + ZT 



(8) 



with w = maxj./ E[Zj;]^^ and z = ey^(maxjy E[Zj,]). 
By expanding Z^^, we calculate E[Z|J and obtain: 

Thus, w < -^=B^B^ and z = ^/2eBuBy. Wc finally obtain the desired result by replacing the expressions of w 
and z in Eq.®. 



Appendix B: Proof of Theorem |4.1 

The proof is composed of 2 steps. 

1. We first perform a simple calculation that will be needed in the second part of the proof. 



Let {Aj}'-^i denote an estimated set of supports having k incorrect atoms: 



J2\A*nA,\ = k. 

We have the equality: 

J J J J 

E E % • ^ " E E • = E E • <^ - E E % • 



j = l (^gA* 



3 = 1 



ipeA*nAj 



as correct atoms in the estimated supports cancel out. Condition in Eq.Q, together with the positivity 
assumption in Eq.([5]) imply that: 



J J 

E E yj" • ~ E E yr'p> vkmy 

j=i -peA] j=i 



(9) 
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2. For a fixed vector of transformations T, we define the set of supports {Aj in tlie following way: 

• as the set of S atoms maximizing dT[i] = X)j=i ■ AjTj{ipi) . 

• For j>2, Aj = Tj{Al) 

By definition, the support Ai is composed of the S atoms maximizing dp. Hence, Ai = A^. 

We say that the estimated supports are h-incorrect, if the total number of incorrectly estimated atoms (in 
all signals combined) is at least equal to h: 

J2\^*nA}\>h, 

We consider the event { JT estimates /i— incorrect supports } and we write the following equalities and 
inclusions on the events: 



{JT estimates /i-incorrect supports} = < (Ai,T) = argmax *^(Ai,r) verifies > A* n A^ 



> h 



= \ t = argmax *s(Af, T) verifies X! I^J* ^ 



Ter 



> h 



C 3r e r verifying *,(A^T*) < ^,{Aj,T) and ^ I A* n A; 



> h 



C U ^vl>,(At,T*)<*«(Af,T),^|A*nAj 



Ter 



>h' 



(10) 



Thus, 



' (JT estimates /i-incorrect supports) < P |J < *s(A^ , T*) < *s(Af , T), ^ I A* n Aj 



> h 



< ^ P U,(At,T*) < ^',(Af ,T),^ |a* n A J 



TeT 



> h 



SJ 



< E E P U.(At,T*)<vI/,(Af,T),^|A*nAj|=fc 



TeTk^lh] 



A*nAj 



From the first part of the proof, wc know that if ^ 
where \E'j, is defined in Section 3. Hence, the following inclusion holds: 



k, then, '^y{Al,T*)-^/y{Al,T) > r^kniy, 



vl/,(AI,T*) < M',(Af ,r),^ |a* n Ajl = fc 



c ,(Aj,r*) < ^y{Ai^T*) - ^1 y |vi/,(Af ,r) > ^y{AlX) + ^} • 
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Hence, 



^s{Al,T*) < ^,(Af , T), Wj n A, 



= k 



(11) 



< 



rjkmy 



rjkmy 



< P 3^ e A*i , X • T* (^) - • A,T* {if) > 



■qkrUy 
25 

rjkmy 
2S 



(12) 
(13) 



We rewrite probability in Eq.(|12| and apply Lemma 4.1 



U ^Xy,.7;*(vp)-s,-A,r,*((p)> 



25 



Similarly for Eq.(113|, 

J 



rjkmy 
2S 



< 2K exp 



Erf] 



4Ci52j + 2C2?7fc5M| J 



> 



rjkmy 
2SJ 



< 



> 



rjkmy 
2SJ 



<2K I- ^^^^^ "^y 

- '''''' I "4Ci52j + 2C72r7fc5M| 



Thus, by combining both results, we obtain: 



SJ 



' ( JT estimates /i-incorrect supports) < ^ ^ 4if exp j — 

Terfe=r/i] \ 



ACiS^J + 2C2i]kSMl 



^^^^^'^'^-P h 4Cr5V + W^ M| 



Let ft, = aSJ, then, we have: 

P (JT estimates (aS J)-incorrect supports) 



/ Number of errors 



SJ 

<l-a) 

< 4:SJK\r\ exp ( 



> a 



MJa^rj'^ ml " 
4Ci + 2C2rja M| 

2 \ 



/ m 

< 4SJK\T\ exp l-CMJa'^rj^j^ 



with C = (4Ci + 26*2 



-1 _ / 32e 



probability on the complementary event. 



§is 



4e\/2) ■ The result of Theorem 4.1 is finally obtained by taking the 
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